from unittest import TestCase
from rcsj_model_analysis import simulating_data as sd

import numpy as np
from scipy import stats
from scipy import constants
from scipy.integrate import quad, simps, cumtrapz
import pymc as pm
import matplotlib.pyplot as plt
import seaborn as sns
import pytensor.tensor as pt
import arviz as az

h_bar = constants.hbar
k_B = constants.Boltzmann
elementary_charge = constants.elementary_charge

        
class TestHazardProperties(TestCase):
    def test_hazard(self):
        I = np.linspace(0, 0.08, 100)
        critical_current = 0.08
        f_j0 = 10**9
        Temp = 30
        
        current = I
        eta = current/critical_current
        freq_j = f_j0*(1-eta**2)**(1/4)
        E_J = h_bar*critical_current/(2*elementary_charge)
        Delta_U = 2*E_J*(np.sqrt(1-eta**2)-eta*np.arccos(eta))
        test_array1 = freq_j*np.exp(-1*Delta_U/(Temp*k_B))
        
        array1 = sd.hazard(I, critical_current, f_j0, Temp)
        
        self.assertEqual(array1.shape,test_array1.shape)
